######################
#  Replication code for 'Mediating the Electoral Connection', forthcoming in the JOP
#  John Henderson and John Brooks
#  12/7/2015    
######################    

# figureII.R
#  :: runs bootstrap sample inference and produces figure II which displays bootstrap IV estimates as evidence of 
#   rain being a strong instrument

rm(list=ls())
setwd('~/Dropbox/rainReplication')

# loop over data and fixed effects specifications 
       
fes.type = 3  
non.missings = 4

source('prelimMain.R')

#detach(covs)
# ANALYSIS             

  
main_iv_fe_dose_full=main_iv_fe_vote_full=main_iv_fe_dose_fes=main_iv_fe_vote_fes=array(NA,nrow(rain_data))
vts=lm(vote~as.factor(fe_id_num),subset=full,data=covs)
dse=lm(dose~as.factor(fe_id_num),subset=full,data=covs)
          
main_iv_fe_dose_fes[as.numeric(names(dse$res))]=dse$res
main_iv_fe_vote_fes[as.numeric(names(vts$res))]=vts$res

vts=lm(vote~as.factor(fe_id_num)+dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + redistricted + 
	dose_prv + vote_prv,
	subset=full,data=covs)
dse=lm(dose~as.factor(fe_id_num)+dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + redistricted + 
	dose_prv + vote_prv,
	subset=full,data=covs)
      
main_iv_fe_dose_full[as.numeric(names(dse$res))]=dse$res
main_iv_fe_vote_full[as.numeric(names(vts$res))]=vts$res
       
main_iv1_fe_fes=ivreg(main_iv_fe_vote_fes~  
	d_inc+
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + redistricted + 
	dose_prv + vote_prv + 	
	main_iv_fe_dose_fes,
	~
	d_inc+
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + redistricted + 
	dose_prv + vote_prv +
	rain_day+rain_day_prev,
	subset=full,data=covs) 

main_iv2_fe_fes=ivreg(main_iv_fe_vote_fes~
	d_inc+
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + redistricted + 
	dose_prv + vote_prv +
	main_iv_fe_dose_fes,
	~
	d_inc+
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + redistricted + 
	dose_prv + vote_prv +
	rain_weekend+rain_weekend_prev,
	subset=full,data=covs)
   
main_iv1_fe_full=ivreg(main_iv_fe_vote_full~
	main_iv_fe_dose_full,
	~rain_day+rain_day_prev,
	subset=full,data=covs) 

main_iv2_fe_full=ivreg(main_iv_fe_vote_full~
	main_iv_fe_dose_full,
	~rain_weekend+rain_weekend_prev,
	subset=full,data=covs)
  
# bootstrap using above residualized outcomes ....

main_iv2_full_boot=array(NA,1000)
main_iv1_full_boot=array(NA,1000) 
  
main_iv2_fes_boot=array(NA,1000) 
main_iv1_fes_boot=array(NA,1000)

  
set.seed(1005)
for(i in 1:1000){
	indx=sample(full,replace=T,size=length(full)) 
	
	main_iv1_fe_fes=ivreg(main_iv_fe_vote_fes~  
		d_inc+
		dist_prev + midterm + pres_party + 
		black + construction + educ + 
		minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
		urban + retail + sos + gov + comp_cq + redistricted + 
		dose_prv + vote_prv + 	
		main_iv_fe_dose_fes,
		~
		d_inc+
		dist_prev + midterm + pres_party + 
		black + construction + educ + 
		minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
		urban + retail + sos + gov + comp_cq + redistricted + 
		dose_prv + vote_prv +
		rain_day+rain_day_prev,
		subset=indx,data=covs) 
   
	main_iv2_fe_fes=ivreg(main_iv_fe_vote_fes~
		d_inc+
		dist_prev + midterm + pres_party + 
		black + construction + educ + 
		minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
		urban + retail + sos + gov + comp_cq + redistricted + 
		dose_prv + vote_prv +
		main_iv_fe_dose_fes,
		~
		d_inc+
		dist_prev + midterm + pres_party + 
		black + construction + educ + 
		minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
		urban + retail + sos + gov + comp_cq + redistricted + 
		dose_prv + vote_prv +
		rain_weekend+rain_weekend_prev,
		subset=indx,data=covs)
                         
	main_iv1_fes_boot[i]=main_iv1_fe_fes$coef[which(names(main_iv1_fe_fes$coef)=='main_iv_fe_dose_fes')]
	main_iv2_fes_boot[i]=main_iv2_fe_fes$coef[which(names(main_iv2_fe_fes$coef)=='main_iv_fe_dose_fes')]
	
	main_iv1_fe_full=ivreg(main_iv_fe_vote_full~
		main_iv_fe_dose_full,
		~rain_day+rain_day_prev,
		subset=indx,data=covs) 

	main_iv2_fe_full=ivreg(main_iv_fe_vote_full~
		main_iv_fe_dose_full,
		~rain_weekend+rain_weekend_prev,
		subset=indx,data=covs)
      
	main_iv1_full_boot[i]=main_iv1_fe_full$coef[which(names(main_iv1_fe_full$coef)=='main_iv_fe_dose_full')]
	main_iv2_full_boot[i]=main_iv2_fe_full$coef[which(names(main_iv2_fe_full$coef)=='main_iv_fe_dose_full')]
   
}            
      
save(
main_iv1_fes_boot,main_iv2_fes_boot,
main_iv1_full_boot,main_iv2_full_boot,
file=paste('bootstrap/figureII-',non.missings,'_',fes.type,'.Rdata',sep=''))
#file=paste('ivResults-bootstrapSEs-',non.missings,'_',fes.type,'.Rdata',sep='')) 

rm(
main_iv1_fes_boot,main_iv2_fes_boot,
main_iv1_full_boot,main_iv2_full_boot
)
        
#######################
# produce plots ...    
 
rm(list=ls())
fes.type=3
non.missings=4

load(paste('bootstrap/figureII-',non.missings,'_',fes.type,'.Rdata',sep=''))

m2=mean(main_iv2_full_boot)
s2=sd(main_iv2_full_boot)
m1=mean(main_iv1_full_boot)
s1=sd(main_iv1_full_boot)

  
set.seed(1006)
x2=rnorm(mean=m2,sd=s2,n=500)
x1=rnorm(mean=m1,sd=s1,n=500)
    
pdf('bootstrap/bootstrap_iv2.pdf')          
plot(density(main_iv2_full_boot),xlim=c(-4,1),ylim=c(0,.8),xlab='Bootstrap IV Estimates',main='')
lines(density(bw=.21,x2),xlim=c(-4,1),lty=2)  
text(x=-3.25,y=.6,label=paste('p = ',0.285,sep=''),cex=1.25)
ks.test(x2,main_iv2_full_boot)
dev.off()

pdf('bootstrap/bootstrap_iv1.pdf')                                
plot(density(main_iv1_full_boot),xlim=c(-4,1),ylim=c(0,.8),xlab='Bootstrap IV Estimates',main='')
lines(density(bw=.21,x1),xlim=c(-4,1),lty=2) 
text(x=-3.25,y=.6,label=paste('p = ',0.212,sep=''),cex=1.25)   
ks.test(x1,main_iv1_full_boot)
dev.off()             

# END      